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The  knowledge  of  physical  properties,  such  as  the  thermal  conductivity,  plays  an  important  role  in  the 
management  of  the  heat  transfer  in  fibrous  materials  like  PEFC  gas  diffusion  layers.  Measurement  of 
thermal  conductivity  by  experimental  means  is  not  easy  (due  to  the  anisotropy  and  the  high  porosity), 
therefore  the  availability  of  experimental  data  is  rather  limited.  In  this  paper,  the  numerical  determina¬ 
tion  of  the  effective  thermal  conductivity  of  fibrous  materials  is  investigated  using  a  three-dimensional 
approach.  Two  different  fiber  geometries  were  studied  with  randomly  generated  fiber  structures  with 
overlapping  and  non-overlapping  fibers.  The  corresponding  anisotropic  thermal  conductivities  are  com¬ 
puted  through  the  solution  of  the  energy  transport  equation.  The  results  were  validated  through  a 
comparison  with  existing  experimental  data  and  the  influence  of  different  parameters  such  as  fiber 
orientation,  fiber  diameter  and  binding  material  were  investigated. 

©  2009  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

The  solution  of  the  purely  diffusive  steady  state  form  of  the  heat 
transfer  equation  requires  the  knowledge  of  one  physical  property, 
namely  the  thermal  conductivity.  This  property  can  be  measured 
experimentally  with  different  techniques  [1],  each  depending  on 
the  type  of  material  or  the  temperature  range,  and  a  high  accu¬ 
racy  can  be  achieved  in  most  of  the  cases.  However,  there  are 
cases,  such  as  fibrous  materials,  where  the  thermal  conductivity 
is  difficult  to  measure  with  high  accuracy  because  of  the  stochas¬ 
tic  nature  of  such  materials  and  the  strong  influence  of  thermal 
contact  between  the  fibers  [2].  Many  parameters  can  have  a  strong 
influence  on  its  value  such  as  porosity,  fiber  geometry  (like  diam¬ 
eter  or  shape),  temperature  range,  fiber  arrangement  (isotropic  or 
anisotropic),  contact  between  fibers,  ratio  between  matrix  and  fiber 
thermal  conductivities,  etc. 

In  order  to  overcome  these  experimental  difficulties,  numerous 
methods  have  been  developed  to  predict  the  thermal  conductivity 
such  as  analytical  means,  statistical  methods,  theoretical  models, 
or  numerical  simulations  [3-7].  However  each  of  these  meth¬ 
ods  has  drawbacks,  like  idealized  structures,  empirical  constants 
(adjustable  parameters),  or  the  need  of  heavy  computing  resources. 
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In  this  work,  the  effective  thermal  conductivity  of  fibrous 
materials  has  been  evaluated  by  using  a  numerical  model  that 
minimizes  the  required  computational  resources.  Due  to  the  three- 
dimensional  character  of  these  simulations,  several  parameters 
can  be  taken  into  account  and  the  full  tensor  of  the  thermal 
conductivity  can  be  determined.  The  method  was  applied  to  gas 
diffusion  layers  of  PEM  fuel  cells  where  the  knowledge  of  the  ther¬ 
mal  conductivity  is  essential  for  the  thermal  management  of  such 
devices. 

2.  Numerical  method 

The  computation  of  the  effective  thermal  conductivity  of  fibrous 
materials  requires  the  solution  of  the  steady  state,  purely  dif¬ 
fusive  three-dimensional  heat  transfer  equation.  In  the  present 
approach,  the  convection  and  radiation  transport,  as  well  as  ther¬ 
mal  contact  resistance  and  phase  changes  were  not  taken  into 
account.  The  solution  can  be  achieved  by  using  several  numeri¬ 
cal  schemes,  such  as  finite  differences  or  finite  elements  [3-5],  by 
analytic  methods  [6],  or  by  fractal  method  [7].  However,  in  the 
case  of  large  three-dimensional  geometries  generated  randomly 
or  based  on  large  data  sets  from  digital  imaging  (like  computer 
tomography),  partial  differential  equation  solvers  are  not  efficient. 
The  description  of  these  complex  geometries  combined  with  the 
large  difference  in  thermal  conductivity  of  the  matrix  and  fibers 
demand  a  very  fine  mesh,  thereby  requiring  large  amounts  of  mem¬ 
ory  and  CPU  time.  Lately,  lattice  Boltzmann  methods  [8-11]  have 
proved  to  overcome  these  limitations  and  can  be  very  efficient 
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Nomenclature 

nx,  ny,  nz  grid  size 
h  size  of  a  voxel  (pan) 

Vf  fiber  volume  fraction 

I<  average  value  of  the  effective  thermal  conductivity 

(WmK-1) 

Kc  effective  thermal  conductivity  of  PCM44  with  car¬ 
bon  fibers 

Kpcm  effective  thermal  conductivity  of  PCM44 


when  working  on  large  geometries.  Weigmann  and  Zemitis  [12] 
used  a  different  approach,  solving  the  energy  equation  by  har¬ 
monic  averaging  and  explicitly  introducing  discontinuities  across 
the  material  interfaces  as  additional  variables.  These  extra  vari¬ 
ables  can  be  determined  through  the  continuity  of  heat  fluxes  at 
the  interfaces.  Fast  Fourier  transform  and  BiCGStab  methods  are 
then  used  to  solve  the  Schur-complement  formulation  for  the  new 
variables.  The  heat  transfer  equation  in  the  purely  diffusion  form 
can  be  written  as 

V.(/CVT)=/  in  Q  (1) 


Let  us  now  consider  a  three-dimensional  parallelepipedic 
domain  of  rq  x  n2  x  n3  cubic  voxels  of  size  #i3.  With  t  the  discrete 
approximation  of  T,  the  harmonic  averaging  discretization  of  Eq. 
(1 )  on  this  grid  is  then 

1  ( (  1  |  1  A  ti+l,j,k-ti,j,k  _  (\  1  A 

h  y\2I<i+ Uk  ^  2KU j,  k  )  h  \2  KiJtk  ^  2  Kt_Uk  ) 

i,j,k  ~  h-1  ,j,k  f  1  1  \  f/,j+l,/c  —  U ,j,k 

x  S  +  h 

V  2KU.k  2 Kij_hk )  h  \  2KiJ,M  2KUj,k ) 

khk+t  -  tjj.fc  _  (  1  \  \  1  tjj'k  -  tjjjc-t  ) 

'  h  \2KiJtk  2Kuhk_, )  h  J 

Ki,j,k  ( ~  ( Ki+1,j,k  ~  Ki,j,k  Ki,j,k  ~  Ki-l,j,k\ 

h  \  Ki+hlk + if,- j, +  i<lhk + if Jik ; 

«  f  KiJ+l,k  -  KiJ,k  KiJ,k  - 
21  \KiJ+i ,k  +  I<ij,k  Kij.k  +  Kij^k  ) 


where  T  is  the  temperature,  I<  is  the  local  conductivity,  / 
represents  the  heat  sources  or  sinks  and  £2  is  the  domain  consid¬ 
ered. 

Although  the  thermal  conductivities  of  the  materials  consti¬ 
tuting  the  fibrous  media  are  considered  in  general  as  isotropic, 
the  effective  thermal  conductivity  of  the  fibrous  material  will  be 
anisotropic  due  to  the  anisotropy  of  the  geometry  and  it  is  fully 
described  by  the  effective  thermal  conductivity  tensor: 
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The  determination  of  the  different  components  of  this  tensor  can 
be  performed  using  the  homogenization  theory.  Three  different 
singular  periodic  boundary  value  problems  have  to  be  solved: 

V  .  (VT„  +  e„))  =  0,  x  e  Q,  (2) 


,  «  I  Ki,j,k+ 1  “  Ki,j,k  ,  Ki,j,k  ~  Ki,j,k- 1  \  \ 

^  Uij.fc+l  +Kfk  KiJtk  +  KiJik_ ,  J  J 

with  /  =  (1 ,  2,  3)  representing  the  direction  of  the  heat  flux.  Since 
the  boundary  conditions  are  supposed  to  be  periodic,  t0  and  tn.+ 1 
can  be  replaced,  respectively  by  tn.  and  t\ . 

For  simplification,  the  derivation  via  jump  conditions  will  be 
described  in  the  following  part  only  for  i  direction.  The  jump  in  the 
heat  flux  [/C(9T/9x i )]  j  k  is  discretized  as 

ts  Tmj,k  ~  Ti+(l/2),jk  T/r  Ti+(l/2)j,k~Ti,j,k 
h/2  iJ’k  #1/2 

=  -8u(Ki+^,k-Kij,k)  +  0(h) 

and  a  first  order  discretization  of  the  condition  of  jump  in  the  tem¬ 
perature  gradient  follows  as 


Tn{x  +  idiei  +  jd2e2  +  kd3e3 )  =  T„(x) 

withn  =  (l,2,3),  ^  =  (0,6#!)  x  (0 ,d2)  x  (0,d3)and  en  is  the  unit  vector 
in  the  respective  direction.  Having  solved  (2),  the  components  of 
the  averaged  conductivity  tensor  can  be  computed  by  integrating 
the  inner  product  of  the  evaluation  direction  i  against  the  flux  and 
vector  in  direction  j: 

K*-^//"K(i)(v,j+5i»‘K  131 

The  difficulty  in  this  approach  consists  in  solving  (2):  the  heat  flux 
I<\/Tn  is  continuous  only  in  directions  perpendicular  to  en  and  has 
a  discontinuity  (or  jump)  in  the  direction  en  that  is  proportional  to 
the  difference  in  conductivities  between  the  matrix  and  the  fibers. 
These  “jump”  conditions  can  be  also  written  as 

[Tn]x=x*  —  0 


K 


din 

dXj 


=  SinlI<k=x*(i  =  1,2,3) 


where  [  ]*=**  is  the  jump  across  the  interface  at  that  position,  dxj 
is  the  directional  derivative  in  the  ith  direction  and  Sin  is  the  Kro- 
necker  symbol. 


8T 
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Using  Taylor  expansions  in  the  presence  of  jumps,  one  finds 

v  i  (Ti+1’J’k  ~  (h/2)[dXl  -  TUj,k 


UfK  \  h 
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After  a  division  by  I<qj0  the  discretization  of  (2)  with  explicit  jumps 
can  be  summarized  as  the  discretization  of  the  partial  differential 


J-iJd-(TIJ-k_1-(fl/2)[a,3T]..N1/2))\ 

h  J 
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equation: 

ti+1  J,k  +  ti-hj,k  I-9*1  0  i+(l/2),J,fc  +  l9x'  i-(l/2)j,k 

h2  2h 
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P  2h 

tij./c+l  +  ti,j,k- 1  LM  ij,k+(l/2)  +  L^J  y, Ml/2) 

+  h2  2h 


6ti 


i,j,k 


/l2 


=  0 


(4) 


and  of  the  jump  conditions  (written  here  for  the  i  direction  for 
simplification): 


2Ki+l,j,k  ~  Kj,j,k  h+1  ,j,k  ~  tjj,k 

^i+ 1  ,j,k  +  ^ij,  k  h 


+  f  ]/_)_( i /2),j,k 


=  -2  8U 


Ki+l,j,k  ~  Ki,j,k 
Ki+1J,k  +Ki,j,k 


Fig.  1.  Fluctuation  of  the  effective  thermal  conductivity  around  the  average  value  I< 
for  three  different  grid  sizes  (n  =  64, 128,  256). 


2 ^i,j,k  Ki-l,j,k  ti,j,k  h-l,j,k 

Ki,j,k  +  Ki-ij,k  h 


=  -2*U 


Kj,j,k  -  Ki-l,j,k 
Ki,j,k  +  Ki-l,j,k 


The  structure  of  (4)  can  also  be  rewritten  in  the  following  way: 


A  \)/ 

V 

V 

D  I 

j 

F 

where  T  is  the  vector  of  temperatures  and  J  is  the  vector  of  nontrivial 
jumps.  The  Poisson  problems  associated  to  this  system  can  then 
be  solved  efficiently  using  FFT,  and  J  is  solved  iteratively  using  a 
BiCGStab  solver. 

Finally,  the  tensor  of  the  effective  thermal  conductivity  can  be 
obtained  by  approximating  (3): 

(h+l,m,p  —  h-l,m,p  l_^Xl  /+( 1/2), m,p  l_^Xl  l-(l/2),m,p  (  ^  ^  ^3 

2 h  4  +°ij )  n 


This  method  is  implemented  in  an  integrated  simulation  environ¬ 
ment  named  GeoDict  [13].  GeoDict  was  used  in  the  current  study 
to  generate  the  structures  and  compute  the  effective  thermal  con¬ 
ductivity. 


3.  Results  and  discussion 


3.1.  Influence  on  numerical  parameters  on  the  results 

Three-dimensional  structures  were  generated  by  specifying 
geometric  parameters  such  as  fiber  diameter  and  length,  fiber  vol¬ 
ume  fraction  and  the  size  of  the  grid  (%xnyx  nz).  As  expected,  the 
generation  of  random  structures  cannot  result  in  a  single  unique 
thermal  conductivity  tensor.  The  values  fluctuate  around  an  aver¬ 
age,  and  the  fluctuation  is  highly  dependent  on  the  size  of  the  grid:  a 
higher  number  of  fibers  will  lead  to  a  smaller  fluctuation.  This  result 
is  illustrated  in  Fig.  1  where  more  then  500  simulations  were  con¬ 
ducted  for  three  different  grid  sizes  (nx  =  ny  =  nz  =  64, 128  and  256). 
For  each  n,  an  average  over  the  500  simulations  was  performed  and 
the  dispersion  of  the  conductivities  around  this  average  was  plot¬ 
ted.  For  the  smaller  grid  (i.e.  n  =  64)  more  than  30%  of  the  structures 
show  a  deviation  (larger  than  9%)  from  the  overall  average  value. 


This  number  decreases  to  10%  and  6%,  respectively  for  the  larger 
grids  (n  =  1 28  and  256,  respectively).  The  influence  of  the  voxel  size 
h  has  also  been  studied,  and  the  value  of  0.5  p,m  was  chosen  in 
the  rest  of  the  study:  a  bigger  value  has  a  negative  effect  on  the 
precision  of  the  results  (the  fibers  are  not  described  with  enough 
precision)  and  a  smaller  value  leads  to  too  small  structures. 

Similarly,  for  a  fixed  grid  size,  a  change  in  the  fiber  volume  frac¬ 
tion  will  have  an  influence:  the  higher  the  fraction,  the  more  fibers 
are  included  in  the  simulation  and  the  smaller  is  the  fluctuation 
(Fig.  2).  Flowever,  the  volume  fraction  size  is  intrinsic  to  the  fibrous 
material  and  cannot  be  varied  to  ensure  more  precise  numerical 
results. 

In  the  following  results,  three-dimensional  structures  were 
generated  using  a  grid  size  of  n  =  256  and  each  particular  set  of 
parameters  was  run  200  times  in  order  to  ensure  high  accuracy  of 
results. 

3.2.  Comparisons  and  validation  with  experimental  data 

Frusteri  et  al.  [14]  studied  experimentally  the  influence  of 
carbon  fibers  on  the  thermal  conductivity  enhancement  of  an 
inorganic  phase  change  material  (PCM44).  Carbon  fibers  of  6  p,m 
diameter  were  added  to  the  phase  change  material  and  the  effective 
thermal  conductivity  was  measured  for  different  volume  fraction 
loadings  and  different  carbon  fiber  lengths  for  a  range  of  aspect 
ratios  (length/diameter)  between  50  and  1 000.  In  the  current  study, 


Fig.  2.  Fluctuation  of  the  effective  thermal  conductivity  around  the  average  value  I< 
for  three  different  fiber  volume  fractions  (Vf  =  0.15,  0.20,  0.25;  n  =  128). 
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Fig.  3.  Evolution  of  the  effective  thermal  conductivity  of  PCM  material  with  carbon 
fiber  loading. 


the  influence  of  the  fiber  length  was  not  investigated  since  it  would 
have  required  a  very  large  grid  (n>  500)  due  to  the  large  range  of 
aspect  ratios.  The  comparison  for  different  volume  fraction  loadings 
between  experimental  results  and  the  current  numerical  results 
is  illustrated  in  Fig.  3.  A  good  agreement  can  be  observed  for  all 
volume  fraction  loadings  for  this  case. 

Khandelwal  and  Mench  [15]  measured  the  through-plane  con¬ 
ductivities  of  Toray  carbon  papers  without  the  addition  of  PTFE  and 
compared  them  with  the  manufacturer  specifications  [16].  These 
results  and  the  numerical  predictions  of  the  current  work  are  sum¬ 
marized  in  Table  1  and  they  appear  to  be  in  good  agreement.  This 
comparison  differs  from  the  previous  one  by  a  higher  volume  frac¬ 
tion  loading  of  carbon  fibers,  and  the  matrix  composition  (air  vs. 
PCM). 

3.3.  Parametric  study 


^  io 
* 

E  9 


0  5  10  15  20  25 

diameter  (|jm) 

Fig.  4.  Evolution  of  the  effective  thermal  conductivity  with  carbon  fiber  diame¬ 
ter  (JCfiber  =  120WmK-1;  I<meAia  =  0.5  WmK-1 ;  porosity  =  0.8— isotropic  distribution 
of  fibers). 


Fig.  5.  Evolution  of  the  effective  thermal  conductivity  with  the  porosity  for  isotropic 
and  anisotropic  fibers  arrangement. 


3.3 A.  Fiber  radius,  porosity  and  fiber  arrangement 

Fig.  4  shows  the  influence  of  the  fiber  diameter  on  the  effective 
thermal  conductivity  at  a  constant  porosity  of  0.8,  with  an  isotropic 
fiber  arrangement.  For  a  diameter  up  to  8  pan,  the  effective  thermal 
conductivity  remains  constant  and  then  decreases  slowly  for  larger 
values.  This  effect  can  be  explained  that  at  a  constant  porosity,  the 
increase  in  diameter  implies  fewer  fibers,  therefore  less  contact 
between  them. 

The  influence  of  the  porosity  for  an  isotropic  and  anisotropic 
arrangement  of  fibers  is  illustrated  in  Fig.  5.  As  expected,  a  lower 
porosity  leads  to  a  higher  value  of  the  effective  thermal  conduc¬ 
tivity  for  both  cases  (isotropic  and  anisotropic).  The  ratio  between 
the  thermal  conductivity  of  the  fiber  over  the  thermal  conductiv¬ 
ity  of  the  matrix  is  over  100:  diminishing  the  porosity  means  more 
material  with  a  higher  conductivity. 

In  the  anisotropic  case,  all  the  fibers  lay  parallel  to  the  x-y  plane 
(Fig.  6).  The  effect  of  the  anisotropy  (Fig.  5)  is  stronger  for  a  high 
porosity,  where  the  ratio  between  the  in-plane  and  the  through- 


Table  1 

Through-plane  thermal  conductivity  of  Toray  carbon  paper  (porosity:  0.78). 


Thermal  conductivity  (W  ml(_1 ) 

Manufacturer  specification  [16] 

1.70 

Khandelwal  experiments  [15] 

1.80  ±0.27 

Current  work  (values  averaged 

1.86 

over  200  simulations) 

"9.790 

Tensor: 

9.820 

1.860 

plane  conductivities  is  over  20,  than  for  a  low  porosity  where  the 
same  ratio  is  less  than  5.  The  in-plane  thermal  conductivity  benefits 
from  the  fibers  length  where  there  is  little  heat  resistance  compared 
to  the  through-plane  conductivity  where  the  matrix  conductivity 
plays  an  important  role. 


Fig.  6.  Anisotropic  arrangement  of  fibers. 
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Fig.  7.  3D  fibrous  structure  with  30wt.%  binder. 


3.4.  Binding  effect 

The  carbon  fibers  constituting  the  gas  diffusion  layers  of  a 
PEM  fuel  cell  are  often  pre-coated  with  polytetrafluoroethylene 
(PTFE),  which  has  a  knock-on  effect  of  changing  the  material  from 
hydrophilic  to  hydrophobic.  PTFE  (commonly  known  as  Teflon) 
has  two  effects:  it  helps  avoiding  water  flooding,  and  it  acts  as 
a  binder  to  maintain  the  integrity  of  the  material.  Manufacturers 
commonly  coat  the  carbon  fibers  with  1 0-30  wt.%  of  PTFE  [  17-20]. 
The  presence  of  this  coating  has  an  effect  on  the  effective  thermal 
conductivity.  Three-dimensional  fiber  structures  were  generated 
with  different  binder  contents  (Fig.  7)  and  their  respective  through- 
plane  effective  thermal  conductivities  were  calculated  and  are 
presented  in  Fig.  8:  the  numerical  results  are  close  to  the  value 
of  1.7  WmK-1  reported  by  the  manufacturer  of  the  Toray  carbon 
paper  (TGP-H-060  without  PTFE).  The  physical  properties  of  the 
different  materials  considered  in  these  simulations  are  presented 
in  Table  2.  The  decrease  of  the  effective  thermal  conductivity  with 
increasing  amount  of  Teflon  can  be  explained  as  follows.  The  PTFE 
acts  as  an  insulator  around  the  carbon  fibers  and  reduces  the  direct 
contact  carbon  to  carbon.  Since  its  thermal  conductivity  is  almost 
500  times  smaller  than  that  one  of  the  carbon  fibers,  it  overall 
reduces  the  effective  thermal  conductivity  of  the  media.  Similar 


Fig.  8.  Evolution  of  the  through-plane  effective  thermal  conductivity  vs.  binder 
loading  (porosity:  0.8). 


Table  2 

Physical  properties  of  carbon  fibers,  air  and  Teflon. 


Thermal  conductivity  (W  mK-1 ) 

Density  (gem-3) 

Carbon  fiber 

130 

1.8 

Air 

0.024 

0.0013 

PTFE 

0.25 

2.2 

behavior  has  also  been  noticed  experimentally  ([1 5,21  ]).  This  effect 
is  well  known  in  the  heat  exchanger  field  where  alternative  com¬ 
posite  coatings  (like  nickel-phosphorous-PTFE)  have  successfully 
replaced  ordinary  PTFE  coatings  [22],  due  to  their  higher  thermal 
conductivity. 

4.  Conclusions 

Three-dimensional  structures  have  been  randomly  generated 
and  the  effective  thermal  conductivity  for  fibrous  materials  has 
been  numerically  determined  under  various  configurations.  Each 
set  of  parameters  was  run  several  times  on  a  large  grid,  ensuring 
higher  accuracy.  Comparison  between  the  results  obtained  from 
this  work  and  experimental  data  obtained  from  literature  and 
from  manufacturers  demonstrated  a  good  agreement.  The  para¬ 
metric  study  demonstrated  that  large  fiber  diameters  reduce  the 
effective  thermal  conductivity  of  the  material,  while  an  increase 
in  the  porosity  leads  to  an  almost  linear  decrease  in  the  ther¬ 
mal  conductivity.  The  anisotropy  effect  (in-plane  vs.  through-plane 
conductivity)  is  more  pronounced  for  large  porosity,  where  the 
contacts  within  fiber  (which  lay  horizontally)  are  reduced.  Finally, 
the  binding  of  the  fibers  with  a  material  of  low  conductivity  (like 
PTFE)  has  a  decreasing  effect  on  the  effective  thermal  conductiv¬ 
ity. 
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